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The study of light-matter interaction has seen a resurgence in recent years, stimulated by highly 
controllable, precise, and modular experiments in cavity quantum electrodynamics (QED) Q]. The 
achievement of strong coupling EHH, where the coupling between a single atom and fundamental 
cavity mode exceeds the decay rates, was a major milestone that opened the doors to a multitude 
of new investigations 0 3?. Here we investigate multimode strong coupling (MMSC) [71 E], where 
the coupling is comparable to the free spectral range (FSR) of the cavity, i.e. the rate at which a 
qubit can absorb a photon from the cavity is comparable to the round trip transit rate of a photon 
in the cavity. We realize, via the circuit QED architecture EDGE], the first experiment accessing the 
MMSC regime, and report remarkably widespread and structured resonance fluorescence, whose ori¬ 
gin extends beyond cavity enhancement of sidebands m Our results capture complex multimode, 
multiphoton processes, and the emergence of ultranarrow linewidths. Beyond the novel phenomena 
presented here, MMSC opens a major new direction in the exploration of light-matter interactions. 


Interest in going beyond strong coupling has focused 
on the ultrastrong coupling limit, where the breakdown 
of the rotating-wave approximation results in excitation 
non-conserving terms [T2] . In contrast, the direction 
which we pursue is the simultaneous strong coupling of 
the qubit to numerous modes, leading to qubit medi¬ 
ated mode-mode interactions and many-body physics not 
present in the single mode problem. Unlike the spin- 
boson problem, where the continuum bosonic modes are 
treated as a bath for the qubit, in MMSC the dynamics 
of individual modes and finite time correlations between 
modes are essential. Furthermore, the study of MMSC 
can not rely on the integrability present in the single 
mode and continuum problems m, thus requiring new 
theoretical ideas for its study. 

The closed system is described by the Hamiltonian: 

H — T ^ ^ ^Eo m ci^ci 17l T % m (cr T cr )(n| n T cim ) 

~ m 

(i) 

where u: a is the qubit frequency, a z is the Pauli operator, 
and m represents cavity mode number. Here, aj^ (cr + ) 
and a m (a~) are mode (qubit) raising and lowering op¬ 
erators. The coupling strength of the qubit to the rath 
harmonic is = go \Jmn + 1, where go is the coupling 
rate to the fundamental cavity mode m 

To achieve MMSC, the coupling between qubit and 
cavity must be comparable to the free spectral range, 
which is made possible through use of a long cavity. Pre¬ 
viously, long cavities have been used for novel comb gen¬ 
eration, with a Kerr nonlinearity provided by the bulk 
medium m- In cavity QED, unique multimode plat¬ 
forms have been studied, where degenerate modes are 
coupled via ensembles of atoms [H GET]. In contrast, 
here the nonlinearity is provided by a single qubit. In 
our setup, the coupling to the fundamental mode [IS], 
go , scales with cavity length L as go ~ 1/L. Holding 
constant the qubit frequency, we are interested in the 
nth harmonic of the cavity which is nearly resonant with 
the qubit. The ratio of the coupling for this mode to the 



FIG. 1. Device for achieving MMSC. a, A low funda¬ 
mental frequency, ujq/2^— 92 MHz, is realized by fabricating 
a 0.68 m coplanar waveguide resonator on a 25 mm x 25 
mm sapphire substrate. Capacitive coupling at the input and 
output ports allows coupling of radiation into and out of the 
resonator, at rates ft m /27r ~ 0.5-3 MHz. Red rectangles show 
coupling capacitors. Blue rectangle encloses qubit, b, Sym¬ 
metric interdigitated capacitors couple cavity to transmission 
lines, c, A transmon qubit is capacitively coupled to the cen¬ 
ter pin near the output capacitor, an antinode for all modes. 
SQUID geometry allows for tuning qubit frequency, co a , via 
the flux bias line (shown) or the external magnet [19] . The 
grid-like pattern in ground plane pins flux vortices. Scale bars 
denote (a) 10 mm, (b) 100 jam, and (c) 100 jam. 

FSR then scales as g n /FSR ~ \[L. As there are many 
modes simultaneously interacting with the qubit, many 
modes can be coupled to each other during the typical 
lifetime of an excitation in the system, with the effective 
qubit mediated interaction rate between a pair of modes 
being (. n-m)FSR i lading t° a delocalization of the exci¬ 
tation across many modes. 

Within this rich domain, we first focus our attention 
on resonance fluorescence, or the response of a qubit to 
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FIG. 2. Demonstrating multimode strong coupling, a, Cavity transmission as qubit frequency is tuned with an applied 
flux bias. Vacuum Rabi avoided crossings are visible as the qubit is tuned into resonance with each mode (bright vertical lines 
spaced by 92 MHz). It is apparent that the vacuum Rabi splitting is comparable with the free spectral range (mode spacing), 
thus indicating that MMSC has been achieved. White rectangles are shown in b, c. b-c, The coupling grows with m: 
go V m + 1. A fit to the first manifold of the multimode Jaynes-Cummings yields go ~ 3.75 MHz. 


coherent drive. Strongly driving a qubit in free space re¬ 
sults in emission forming a three peaked structure known 
as the Mollow triplet, composed of a center peak and two 
symmetric sidebands [20]. The sidebands are displaced 
linearly from the central peak by the Rabi amplitude 
[21] Q. This phenomenon, also observed experimentally 
with superconducting qubits [22], is typically explained 
using the dressed state picture approach [23]. Coupling 
of a qubit to a single mode cavity strongly modifies the 
Mollow triplet [24V 26 . For weak coupling, the cavity 
passively filters the fluorescence, while at strong cou¬ 
pling, the sideband width is proportional to the coupling 

nun]- 

We design a very low frequency cavity such that the uj a 
of a standard qubit falls near a high harmonic (^ 50th- 
75th) of the cavity, unlike traditional circuit QED exper¬ 
iments. Additionally, this approach allows us to adapt 
to the typical measurement range (3-8GHz) of existing 
circuit QED setups (see Methods). The superconduct¬ 
ing microwave cavity comprises a 0.68 m transmission 
line with identical capacitors on either end, resulting in 
a fundamental frequency, uoq/2tt = 92MHz (see Fig. 1 and 
Methods section for details). We couple a flux-tunable 
transmon qubit [19! near one end of the resonator, an 
ant i-node for all the modes of the cavity, and achieve 
a qubit-mode coupling strength for the 75th cavity har¬ 
monic exceeding 30 MHz. Coherent drive, at frequency 
is introduced only via the input port (Fig. lb), while 


radiative decay of cavity photons occurs via both input 
and output ports. A unitary transformation maps the 
coherent drive from the cavity to the qubit, yielding a 
qubit driving term Q cos (ujd t)(cr + +cr _ ), with Q the Rabi 
drive amplitude (see Supplement). The mode linewidths, 
ft m /27r, span from 0.5-3 MHz to provide reasonable sep¬ 
aration between modes. The qubit spontaneous decay 
rate, 7/27T, is of the order of 1.6 MHz, due to the multi- 
mode Purcell effect m • 

Measuring transmission while tuning the qubit fre¬ 
quency reveals that the qubit strongly couples to a vast 
number of modes (Fig. 2a) [3]. In fact, we are able to fol¬ 
low the trail of avoided crossings for several GHz as the 
qubit energy tunes through the uniformly dense mode 
structure. The smooth transitions between subsequent 
avoided crossings show that the qubit is perpetually near 
resonance with some mode and as such the dispersive 
approximation is never valid. In Figs. 2b-c, the marked 
difference in the magnitude of the avoided crossings, and 
thus coupling strength, is apparent. We apply a fit to 
find g 0 /27T ~ 3.75 MHz, giving rise to g7^/2n ~ 32 MHz. 

In Fig. 3 we observe fluorescence across more than 50 
modes, when the drive and qubit are resonant with a 
high harmonic of the cavity. Enhanced fluorescence is 
observed at mode m when it is resonant with the Rabi 
sideband. As the displacement of Mollow sidebands is 
proportional to drive amplitude f2, this occurs dl when 
detuning from the rath mode A m = oj rn — uj ( i ~ Q. For 
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FIG. 3. Multimode resonance fluorescence, a, Power spectrum for varying drive power, Pd , when ujd — — <^ 75 , measured 

at more than 50 modes. Mode index n is the relative detuning from the drive mode, n = m — 75. Power is measured in a 
10MHz window around each mode. Increasing drive power unveils fluorescence at farther detuned modes, strong emission from 
multiphoton processes, and complex multi-lobed fluorescence spectra, b, The drive power corresponding to peak fluorescence 
for each far detuned mode shows a quadratic dependence on detuning, suggesting this fluorescence is due to cavity enhanced 
sideband fluorescence nu. Modes near the drive frequency (see inset) do not follow the quadratic fit due to strong multimode 
interaction, c, Detailed plot of n = —30 from (a) shows a multi-lobed structure in fluorescence spectrum. Bare mode frequency 
is denoted with a vertical white dashed line. For d-g, uj d — uo a — ^69- d, We explore the difference between the two lobes seen in 
(c) for several modes (n m —36, —34, —32, —30, —28, from bottom to top with individual curves vertically displaced for clarity) 
by measuring the total fluorescence power in a 6 MHz window while varying drive power. The drive powers corresponding to 
peak fluorescence move oppositely with mode number for the two lobes, e-g, For modes near the drive (n= -1, 1, 2), strong 
multimode interaction further yields a more complex multi-lobed structure. Bare cavity mode frequencies are denoted with 
vertical white dashed lines. 


far detuned modes the drive power needed to reach the 
peak fluorescence at these modes pP eak ^ and hence 
P p d eak » A l (Fig. 3b). 

At fixed drive power, emission is spread over many 
modes, with the simultaneous enhancement at many 
nearby modes, differentiating the multimode fluorescence 
from the single mode case. With varying drive power, far 
detuned modes exhibit a multi-lobed fluorescence struc¬ 
ture (e.g. Fig. 3c). The drive power for the peak of the 
secondary lobe approaches twice that of the first lobe for 
increasing detuning A m (Fig. 3d), originating from mul- 
tiphoton decays. Two photon decay processes in single 
mode cavities, where = 2A m , have previously been 
observed HU [28]. For MMSC, the multiphoton decay 
processes can progress via combinations of many differ¬ 
ent modes. The resonance condition for single and mul¬ 
tiphoton processes is thus shifted away from f 1 = p A m , 
for p G Z (see Supplement). 


At smaller detuning, the drive power for the peak fluo¬ 
rescence deviates from a quadratic dependence on detun¬ 
ing. Multimode effects become increasingly dominant as 
we approach modes closer to the qubit, with the appear¬ 
ance of a different intricate multi-lobed structure. When 
the resonant mode is very weakly driven, there is no in¬ 
coherent scattering to the other modes. However, at a 
higher drive power, fluorescence appears simultaneously 
at multiple nearby modes (three of which are displayed 
in Fig. 3e-g). In addition to detecting bright fluorescence 
at these nearby modes (three on either side of the cen¬ 
ter mode) we also observe incoherent scattering at the 
center mode itself. Furthermore, there is a change in 
emission frequency at each of these modes, approaching 
the bare cavity mode from above with increasing drive 
power. This shift is qualitatively the same for modes on 
either side of the drive, and cannot be attributed to a 
simple AC Stark shift. 
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FIG. 4. Ultranarrowing of fluorescence. a, Power 
spectrum (n — 1) for varying drive frequency (wd/ 27t), at 
fixed drive power, reveals linewidth narrowing. Minimum 
linewidth, 65kHz, corresponds to (ujd — cj69)/27r = 2 MHz. 
b, Normalized spectra along dashed lines in (a) are plotted in 
blue, vertically displaced. For comparison, low power trans¬ 
mission is plotted in solid gray and the bare mode frequency 
is denoted with a vertical gray line, c, Emitted power (over ~ 
3 MHz bandwidth) vs. Ud/ 27r shows simultaneous brighten¬ 
ing of n — =hl modes, d, Fluorescence linewidths vs. Ud/^Ti. 
Linewidths vary concurrently for n = =bl and are inversely 
proportional to emitted power. For comparison, dashed lines 
show k/2tt. 

As drive power is further increased, the simultaneous 
fluorescence at the nearby modes fades away. At an even 
higher power, the emission at these modes reappears, but 
at the bare cavity mode frequencies. This power coincides 
with that at which the far detuned modes fluoresce, sug¬ 
gesting the emergence of qubit mediated massively mul¬ 
timode correlation. 

A remarkable observation of our experiment is the ap¬ 
pearance of ultranarrow linewidths in emission. The 
fluorescence at far detuned modes, captured in Fig. 3, 
have linewidths that are roughly equal to that of the 
bare modes. For modes nearest the qubit the fluores¬ 
cence linewidth is smaller than e.g. for the nearest 
mode above being respectively 650 kHz vs. 1.1 MHz. 
Detuning the drive from the center mode unveils here a 
linewidth narrowing by over an order of magnitude to 65 
kHz (Fig. 4a). The narrowest fluorescence occurs when 
the drive is detuned by 2 MHz from the resonant mode. 
This unexpected behavior is robust to qubit frequency; 
the narrowed fluorescence produced with the 2 MHz de¬ 
tuned drive remains even as the qubit frequency is tuned 
over 150 MHz away from the cavity mode (see Supple¬ 
ment). 

We attribute this narrowing to the spontaneous gen¬ 
eration of coherence arising from the coupling of many 
dressed states through a common vacuum. Such behavior 


has been predicted and observed in atoms with multiple 
closely spaced bare atomic levels coupled to vacuum or 
a single mode of a cavity [29H3T] . where the theoretical 
analysis is aided by the simpler dressed state structure. 
In our setup, unlike the single mode case, the number 
of states in a single excitation manifold is not a constant 
but rather grows as + C^m^-i 2 ) > where M is the 

number of modes and N is the excitation manifold. The 
rapid growth of the number of states per manifold hinders 
a simplified analysis via the dressed state picture, but 
makes possible the near resonance of many level spacings 
(Bohr frequencies), leading to a collective enhancement of 
the coherences from their mutual couplings. An external 
coherent drive allows us to dynamically access regions of 
the dressed states with large generated coherences, as ev¬ 
ident in the drive dependence of the narrowing in Fig. 4. 

In the region where narrowing is observed, the 
linewidth and steady state photon number are observed 
to be inversely proportional, as has been previously sug¬ 
gested for the case of a single mode [32] . For our multi- 
mode cavity, the narrowing is in fact observed simulta¬ 
neously for multiple modes at the same drive frequency, 
as depicted in Figs. 4c-d for the two modes nearest the 
drive. While the modes directly neighboring the cen¬ 
ter mode exhibit the sharpest linewidths of 65 kHz, the 
next nearest modes also display a narrowed fluorescence, 
^ 300 kHz. The multi-lobed power dependent structure 
of fluorescence also appears at the optimal detuned drive 
frequency, however, the linewidth of the second lobe is 
of order /c m , for the nearby modes. For more distant 
modes, while fluorescence follows a similar depen¬ 
dence, there is no narrowing apparent. Finally, the drive 
dependent narrowing is not unique to these described ex¬ 
perimental parameters. When we translate the drive and 
qubit to another set of modes, we observe equally narrow 
linewidth fluorescence which indicates that spontaneous 
generation of coherence is generic in this system. 

In this letter, we realized a new regime of circuit QED, 
multimode strong coupling. We demonstrated the cre¬ 
ation of many-body steady states with long lived quan¬ 
tum coherence, as evidenced by multimode resonance 
fluorescence with ultranarrow linewidths. Dynamically 
generated quantum correlations between modes in this 
driven dissipative system is expected to lead to qualita¬ 
tively new physics [33]. While this new complex regime 
poses many theoretical challenges, precise circuit QED 
experiments can guide this exploration. We hope that 
access to multimode strong coupling will lead to many 
new advancements in circuit QED and our understand¬ 
ing of light-matter interactions. 

Methods 

The MMSC device was made using a combination of 
standard optical and electron-beam lithography tech¬ 
niques on a 25 mm x 25 mm sapphire wafer. The cavity 
and larger features were etched from a 200 nm thin 
film of sputtered Niobium, while the qubit was made 
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with evaporated Aluminum using the Dolan Bridge 
technique. The cavity has a center pin to gap ratio 
of 10:4.186. The transmon qubit [19] is capacitively 
coupled near the end of the cavity, at an antinode. The 
qubit consists of two large islands and a pair of parallel 
junctions in a SQUID-loop geometry. While there is an 
on-chip flux bias line, a superconducting magnet was 
installed around the device to tune the qubit frequency 
without detrimental thermal heating. The device was 
extensively wire-bonded to the circuit board to ensure 
proper connection of all ground planes and reduce the 
effect of spurious modes. The device was cooled to 
11 mK in a dilution refrigerator, and measured using 
standard low-temperature microwave circuit elements 
including an IQ mixer, isolators, and a HEMT amplifier. 
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I. EXPERIMENTAL DETAILS 
A. Fridge Diagram 

The experimental measurement setup is shown in Fig. 

SI. 


B. Standard Spectroscopy 

Two-tone spectroscopy or dispersive readout is a stan¬ 
dard measurement technique which uses the dispersive 
shift of a mode to precisely determine the frequency of 
the qubit |Tj . In our massively multimode strong coupling 
experiment, such spectroscopy shows that the mode fre¬ 
quency is affected not only by the qubit state, but also by 
populating the modes nearby to the qubit. This is yet an¬ 
other manifestation of mode-mode interactions mediated 
via the qubit. The magnitude of transmission, a proxy 
for the dispersive shift, is indicated by the colorbar in 
Fig. S2. As the qubit frequency approaches a mode, 
the strength of interaction between that mode and the 
measurement mode steadily increases. 

This interaction is so strong that the measurement 
mode shifts substantially even when the qubit is over 200 
MHz away. Thus the frequency of the measurement mode 
needed to be recalibrated for each flux voltage in Fig. 
S2. Furthermore, to track the qubit in this frequency 
range with sufficient contrast, the designated measure¬ 
ment mode needed to be changed such that it was far¬ 
ther and farther away from the qubit. For each change 
of measurement mode, there is a horizontal black line in 
Fig. S2. 
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C. Ultranarrow Linewidth Versus Qubit and Drive 
Frequencies 

We find that ultranarrow fluorescence is very robust to 
the flux bias voltage of the transmon qubit as shown in 
Fig. S3, where we display the linewidth as a function of 
the flux bias voltage applied to the qubit. This observed 
robustness will be discussed in the following theoretical 
analysis. Also displayed is the integrated fluorescence 
power against the same bias. The power and linewidth 
are seen to be inversely related. The drive frequency and 
power here are both fixed to that which results in the 
minimum linewidth in Fig. 4 of the main text. 


II. THEORETICAL DISCUSSION 

A. Effective Hamiltonian and Master Equation 

In order to better interpret the intriguing dynamical 
phenomena observed in our experiment, here we provide 
a simplified model to illustrate key aspects of the physics. 
Even though the model is far from a complete represen¬ 
tation of the experimental system, we argue that qualita¬ 
tively similar physics emerges in the massively multimode 
system. 

This driven dissipative system can be described by the 
following master equation (ft = 1 in our units) 


= -i [H, p] + \ ( 2cr P° + - cr+cr P - pcr + cr~) ( 2a ™P a m - a f m a mP - pa f m a m ) , (2) 

rri 

H = H 0 +ip (4 - a r ) , (3) 

Ho = ^2 A m a f m a m + A a o z + g m (4™°"” + h.c.) , ( 4 ) 

rri 


with Km the photon loss rate for mode m, and 7 = 1/Ti 
the qubit spontaneous decay rate. The above Hamilto¬ 
nian is written in the rotating frame of the drive, and 
the rotating wave approximation (RWA) has been used. 
A m — — Ud is the detuning between the rath mode 

and the drive. Likewise, A a = uj a — ujd is the detuning 
between the qubit and the drive. RWA is a good ap¬ 
proximation if we focus on high harmonics 70 in the 
experiment) of this multimode system. 77 is the ampli¬ 
tude of the external drive, and the drive is assumed to be 
resonant with the rth mode. Steady state photon num¬ 
bers of other modes in the uncoupled system (g m = 0 ), 
are suppressed relative to the resonant term, by multi¬ 
ples of (FSR/k ) 2 , and we are justified in keeping only 
the resonant drive. Applying an active drive and dissipa¬ 
tion dependent unitary displacement to the state of the 
system, via p = pD, with 

D = exp (£(4 — a r)) , , (5) 


transfers the drive in the master equation onto the qubit. 
The Hamiltonian appearing in the master equation is now 
H = Hq + Q(cr + + a - )/ 2 , where Cl = Sg r g/K r is the Rabi 
amplitude. We could in fact directly implement such a 
drive in our experiment by pumping the flux bias line 
to qubit. For the sake of simplicity we start from this 
form of the Hamiltonian (externally driving the qubit) in 
what follows, and study a restricted two mode system, 
also taking the drive to be resonant with the qubit. We 
disregard the mode resonant with the qubit, focusing on 
the more interesting sideband phenomena. 


The master equation for the simplified two mode model 
is 


dp_ 

dt 


-i [H, p} + ir (2 a pa + - a + a p ■ 


pcr+cj-) + ^ ~y ( 2 

m= 1,2 


CL rn p0^ rn (l^nCLffiP pOj^dfYi^j , 




+ h.c.^ + g 2 {a\(J + h.c.^ 


(6) 

( 7 ) 


Without cavity modes, the qubit dynamics result in the a special choice of parameters, can be solved analytically 
famous Mollow triplet. While the single mode case, for (Cl = A) |2j, this is in general not possible for a general 






multimode system. We thus resort to numerical simu¬ 
lation based on the above equation ^ in the following 
discussion. 

To motivate the numerical exploration of the two mode 
master equation, we first analyze some properties of this 


Hamiltonian to gain intuition. Introducing a new ba¬ 
sis for the qubit by defining |0) = 4| (|0) — 11)), 11) = 

^ (|0) + |1)), yields the Hamiltonian 


SI 


H = o^ + E A 


t I fh 
a ai + — 


i=1 


(a z - a + + a ) a\ + h.c. 


Applying an atomic-state-dependent mode displacement 
operator U = Y\*i=i ex P ( 2^7 ( a i ~ a \ )°z) to remove o z 


from the interaction, we have 


( 8 ) 


rr si ~ 

H = —G z 
2 


E* 


CL - CLi 


9i 


i=1 


(a+VM ~ V-VM) (“i ~ 4) 


(9) 


where 

Vi = exp ^(aj - a i)^j (10) 

= exp (-(gi/ Ai) * 2 /2) ^ (&/ \) m+n f A af m a?. 

z ' m!n! 

m,n 

(ii) 

This form of the Hamiltonian reveals the complicated 
multimode multiphoton interactions of arbitrary order. 


To further simplify the analysis, as in the RWA, we 
truncate the Taylor expansion to remove fast oscillating 
terms. These approximations are only valid in the limit 
of gi/Ai <C 1, and near certain resonance conditions. 
We thus arrive at a finite order multimode multiphoton 
Hamiltonian. For example, let us assume that Ai = 
—A 2 = A,#i = #2 = g, Ki = k>2 = ft, and that the qubit 
is resonant with the drive frequency. In this special case, 
the resonance conditions are simply SI « mA, for m G 
(1,2). The effective Hamiltonians for the two different 
Rabi rates (Q = A and Q = 2A) become 


H^ f A) = Aajai - A a\a 2 + ^-g z + |[d + ai - d + a]J> + h.c.] , 

2 

H^t- 2 /s) _ _ ^ a t fl2 _|_ cr+ a t + h.c.] . 

2A 


rn = 1 

(12) 

m = 2 

(13) 


Note that two photon processes in the second Hamilto¬ 
nian (13) above are supressed relative to the single pho¬ 
ton ones in ( p"2| ) by gj A. Terms like d+aj arise because 
the second mode carries negative frequency in the rotat¬ 
ing frame. The second Hamiltonian is the basis for two 
photon lasing of driven atoms in a cavity [3] . In our case 
where only a single qubit is present, keeping high-order 
quantum correlations beyond semiclassics is essential in 
explaining relevant coherence properties of the fluores¬ 
cence light. The full master equation could be derived 
by applying these transformations on the dissipative part, 
giving rise to complicated multimode multiphoton dissi¬ 
pative terms. We instead jump directly to the numerical 


results derived from the two mode master equation ©• 


B. Multimode Multiphoton Processes 

We use the Monte Carlo Wavefunction approach to 
numerically simulate the dynamics of this driven dissi¬ 
pative multimode system. In our simulation, the size 
of the Hilbert space needed grows as 2EL N m , with 
7V m the photon cutoff for mode m. The rapid growth 
of the dimension of the Hilbert limits the number of 
modes that can be incorporated into the simulation. Ex¬ 
pectation values of relevant operators are followed and 
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steady state values for the driven system are presented 
below. To compute two time correlation functions such 
as (a\ 2 (t) ai, 2 ( 0 )) using this method, we have to employ 
a two step procedure [5], equivalent to the quantum re¬ 
gression theorem, but at the expense of requiring large 
numbers of walks to produce reliable results. 

To explain the lobe structure and multimode fluores¬ 
cence observed in our experiment, we first analyze two 
special limits of the two mode model. Resonant enhance¬ 
ment by multiphoton processes in a single mode i (via the 
interaction cr + a™ +h.c.), while tuning the qubit drive U, 
leads to peaks in the steady state photon number near 
subharmonics of the Rabi rate |A^|=U/ra,raEZ. This 
is the physical origin of the lobe structure at far detuned 
modes in our massively multimode system, as g/A <C 1 is 
a good approximation for these modes. Fig. S4 presents 
numerical results for the special two mode case where 
the two modes are situated symmetrically around the 
qubit (and the drive which is resonant with it). Flu¬ 
orescence at these two modes are equally and simulta¬ 
neously enhanced. Moreover, the two peaks in photon 
number and qubit polarization cr + deviate from the bare 
resonance conditions |A| = Q/m,m = 1 , 2 . The qubit 
inversion a z can go above 0. The results for a single 
mode model are also shown for comparison. While mul- 
tiphoton resonances are apparent in both cases, the two 
mode system reaches much larger photon numbers in the 
steady state. The difference also appears in the spectral 
linewidth (see Fig. S 6 ). While the single mode fluores¬ 
cence shows no spectral narrowing, the spectral linewidth 
in the two mode model can be much smaller than the nat¬ 
ural linewidth. Beside the simple two photon interaction 
in a single mode (cj~a\a\ ) analyzed above, the interaction 
vertex coupling two different modes (cr _ cqa^,i ^ j ) will 
result in resonant enhancement in modes near harmon¬ 
ics of the Rabi rate |A| = rriQ. This is demonstrated in 
Fig. [5] where another special two mode case (A 2 = 2Ai) 
is investigated. Similarly to the previous case, this effect 
leads to co-enhancement of fluorescence at two modes, 
as evidenced by the brightening of the second mode at 
A ~ 80 MHz. As this is a second-order effect induced by 
the first mode, N 2 is much smaller than Ni. For stronger 
coupling, we would expect the difference to be smaller. 

In our strong coupling multimode system, all multi- 
mode multiphoton terms have to be taken into account, 
which makes the analysis very cumbersome. The compe¬ 
tition between these processes will generally make reso¬ 
nances broader and less prominent. The numerical simu¬ 
lations qualitatively agree with our experimental results. 

C. Multimode Spectral Narrowing 

Ultranarrow resonance fluorescence has been predicted 
for the single mode cavity [4], where coherence between 
dressed states of the driven system arises from coupling 
to common vacuum fluctuations. The resultant spec¬ 
trum linewidth of the mode could be much smaller than 


the natural linewidth of the bare system. In fact, the 
linewidth of the sideband is proportional to the inverse 
of the steady state photon number 1/N. We provide nu¬ 
merical evidence that this is also true in the multimode 
system. It is worth noting that reference [4] only con¬ 
sidered a special value of the parameters, where U = A 
and g ~ 7 . Here we generalize their result to explain our 
observations. 

In Fig. S 6 , for the symmetric two mode system, it is 
shown that spectral narrowing appears for a range of the 
Rabi rate U. Moreover, the linewidth is inversely propor¬ 
tional to the photon number (see also Fig. S4(a)) such 
that the linewidth reaches its minimum as the photon 
number reaches its maximum at U = 106 MHz. The left 
mode also shows the same amount of spectral narrowing 
(data not shown). The spectral narrowing in this simpli¬ 
fied model doesn’t reach that of our experiments (a factor 
of ^ 3.3 in the simulation vs. ^ 17 in the experiment), 
yet it displays qualitively similar features. Moreover, the 
relation Kfi uo ~ 1/N suggests that the linewidth could be 
made much narrower than that of the current simplified 
model, when reaching larger steady state photon num¬ 
bers. For the interacting massively multimode system, 
the proliferation of the dressed states in the high excita¬ 
tion manifolds leads to enhanced generation of coherence 
among the dressed states from a collective interference, 
which relies on the near resonance of many dressed state 
level spacings. Concomitant with the collective suppres¬ 
sion of dissipation for the steady state then is the rise in 
the steady state occupation. Simulations of such a mas¬ 
sively multimode system demands more computational 
resources due to the growth of the Hilbert space. More 
extensive results will be presented elsewhere. 

Spectral narrowing is generally very insensitive to the 
qubit frequency. For finite drive-qubit detuning A a , the 
effective Rabi rate is U = (fig + A 2 ) 1 / 2 , with fi 0 the 
amplitude of the drive. Based on the above analysis, 
spectral narrowing occurs when fi is within the range 
A =b £(fio), where A is again the mode frequency mea¬ 
sured with respect to the drive (i.e. its detuning), and 
£(fio) is generally determined by the number of modes, 
the coupling strength and the dissipation. For relatively 
small detuning A a and large Rabi rate fio, this gives a 
large effective range of qubit frequency S(A a ) ~ -§~S(Q 0 ) 
over which narrowing occurs. It is worth noting that even 
if the drive is not resonant with the qubit and the two 
modes are not symmetric around the drive, spectral nar¬ 
rowing still exists in our simulation. 

For a single mode system, it was argued that spec¬ 
tral narrowing requires that the qubit decay rate 7 be 
comparable to the coupling strength g , i.e. the moder¬ 
ate coupling regime [4]. We therefore study the effect of 
qubit decay rate 7 on steady state properties of the two 
mode system. Fig. S7(a) presents photon number N± 
versus Rabi rate for a different 7 value. From this result, 
we conclude that spectral narrowing survives when going 
from moderate coupling to the strong coupling regime, 
which we attribute to the fact that the eigenstates of the 
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multimode system in high excitation manifolds are more 
densely spaced than in the single mode system. Dynam¬ 


ically generated coherence between these closely spaced 
states then results in spectral narrowing of many modes 
simultaneously. 
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Figure S 1. Schematic representation of measurement setup. Device is mounted within the superconducting magnet, inside a 
dilution refrigerator. Experimental results presented in the paper were obtained from this general setup, with varied digital 
processing. 
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Figure S 2. Mode-mode interaction in dispersive readout. Tracking the change in dispersive shift for measurement mode, at 
cjm, with varying probe frequency, uj p , shows the impact of populating nearby modes. The magnitude of the dispersive shift is 
related to the detuning of a mode from the qubit (at 7.08GHz for <j)b — 0.6). The qubit frequency is tuned using the on-chip 
flux bias line. Fainter curves are attributed to higher order processes. Horizontal black lines indicate a change of ujm (6.469, 
6.653, and 6.747 GHz from top to bottom), chosen to maintain sufficient distance from the qubit. 




Figure S3. Measured linewidth and emitted power versus flux bias. Spectral narrowing appears over a wide range of flux bias, 
corresponding to over 150 MHz of the qubit frequency. The other parameters are the same as in Fig. 4 in the main text. The 
fluorescence power and the linewidth are inversely related. 
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Figure S 4. Steady state photon number, qubit inversion and polarization, first parameter set. Parameters of the two mode 
system are Ai/27r = — A 2 / 27 T = 100 MHz,/^i/27r = — 1 MHz,7/27t = 15 MHz,#i/27r — <72/271- — 15 MHz. Parameters 

of the single mode system are A/27r = 100 MHz,ft/27r = 1 MHz, 7 / 271 - =» 15 MHz,#/27r = 15 MHz. These plots include 24 
different Rabi rates Q with 1,000 quantum walks per Rabi rate in the numerical simulations. 



Figure S5. Steady state photon number, qubit inversion and polarization, second parameter set. Parameters of the two mode 
system are 2 Ai/27t = A 2 / 27 T = 160 MHz, fti/27r = K2/2tv = 1 MHz,7/27r = 15 MHz,gi/27r = g2/2n = 25 MHz. These plots 
include 24 different Rabi rates Q with 1,000 quantum walks per Rabi rate in the numerical simulations. 
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Figure S 6 . Fluorescence spectrum of the right mode in a two mode system, (a) Parameters of the simulation are the same 
as Fig. S4, with Q/ 2tt = 106 MHz, where the photon number reaches its peak. This figure includes 12,000 quantum walks 
in the numerical simulation. The spectrum is generated by Fourier transforming the two time correlation function of the 
mode operator. The fluorescence linewidth is more than 3 times narrower than the natural mode linewidth. (b) Fluorescence 
linewidth versus Rabi rate. The black dashed line indicates the natural linewidth of the mode (1 MHz). 




Q (MHz) Q (MHz) 


Figure S 7. Photon number and linewidth of the right mode for varying drive amplitude. Parameters are the same as Fig. S4 
except that 7 = 5 MHz instead of 15 MHz. The sharpest line does not correspond to the brightest fluorescence. Nevertheless, 
it is a factor of 2 narrower than the natural linewidth of the bare system. The black dashed line indicates the natural linewidth 
of the mode (1 MHz). 










